-- Kuwahara filter
--
-- 
-- Performs the Kuwahara Filter. This filter is an edge-preserving filter.
--
--
-- ( a  a  ab   b  b)
-- ( a  a  ab   b  b)
-- (ac ac abcd bd bd)
-- ( c  c  cd   d  d)
-- ( c  c  cd   d  d)
--
-- In each of the four regions (a, b, c, d), the mean brightness and the variance are calculated. The
-- output value of the center pixel (abcd) in the window is the mean value of that region that has the
-- smallest variance.
--
-- description copied from http://www.incx.nec.co.jp/imap-vision/library/wouter/kuwahara.html
--
-- implemented by Øyvind Kolås <oeyvindk@hig.no> 2004


-- the sampling window is:  width=size*2+1  height=size*2+1

size = 4

edge_duplicate = 1;

-- local function to get the mean value, and
-- variance from the rectangular area specified
function mean_and_variance (x0,y0,x1,y1)
 local variance
 local mean

 local min         = 1.0
 local max         = 0.0
 local accumulated = 0
 local count       = 0
   
 local x, y

 for y=y0,y1 do
   for x=x0,x1 do
     local v = get_value(x,y)

     accumulated   = accumulated + v
     count         = count + 1
      
     if v<min then min = v end
     if v>max then max = v end

    end
  end

  variance = max-min
  mean     = accumulated /count
   
  return mean, variance
end

-- local function to get the mean value, and
-- variance from the rectangular area specified
function rgb_mean_and_variance (x0,y0,x1,y1)
 local variance
 local mean
 local r_mean
 local g_mean
 local b_mean 
 
 local min           = 1.0
 local max           = 0.0
 local accumulated_r = 0
 local accumulated_g = 0
 local accumulated_b = 0
 local count         = 0
   
 local x, y

 for y=y0,y1 do
   for x=x0,x1 do
     local v     = get_value(x,y)
     local r,g,b = get_rgb (x,y)

     accumulated_r   = accumulated_r + r
     accumulated_g   = accumulated_g + g
     accumulated_b   = accumulated_b + b
     count         = count + 1
      
     if v<min then min = v end
     if v>max then max = v end

    end
  end

  variance = max-min
  mean_r   = accumulated_r /count
  mean_g   = accumulated_g /count
  mean_b   = accumulated_b /count
    
  return mean_r, mean_g, mean_b, variance
end

-- return the kuwahara computed value
function kuwahara(x, y, size)
  local best_mean     = 1.0
  local best_variance = 1.0

  local mean, variance
	
  mean, variance = mean_and_variance (x-size, y-size, x, y)	

  if variance < best_variance then
    best_mean     = mean
    best_variance = variance
  end
	
  mean, variance = mean_and_variance (x, y-size, x+size,y)

  if variance < best_variance then
    best_mean     = mean
    best_variance = variance
  end	
	
  mean, variance = mean_and_variance (x, y, x+size, y+size)

  if variance < best_variance then
    best_mean     = mean
    best_variance = variance
  end	
	
  mean, variance = mean_and_variance (x-size, y, x,y+size)

  if variance < best_variance then
    best_mean     = mean
    best_variance = variance
  end	

  return best_mean
end


-- return the kuwahara computed value
function rgb_kuwahara(x, y, size)
  local best_r, best_g, best_b
  local best_variance = 1.0

  local r,g,b, variance
	
  r,g,b, variance = rgb_mean_and_variance (x-size, y-size, x, y)	

  if variance < best_variance then
    best_r, best_g, best_b = r, g, b
    
    best_variance = variance
  end
	
  r,g,b, variance = rgb_mean_and_variance (x, y-size, x+size,y)

  if variance < best_variance then
    best_r, best_g, best_b = r, g, b
    best_variance = variance
  end	
	
  r,g,b, variance = rgb_mean_and_variance (x, y, x+size, y+size)

  if variance < best_variance then
    best_r, best_g, best_b = r, g, b
    best_variance = variance
  end	
	
  r,g,b, variance = rgb_mean_and_variance (x-size, y, x,y+size)

  if variance < best_variance then
    best_r, best_g, best_b = r, g, b
    best_variance = variance
  end	

  return best_r, best_g, best_b
end

function kuwahara(radius)
  for y=0, height-1 do
    for x=0, width-1 do
      r,g,b = rgb_kuwahara (x,y, radius)
      set_rgb (x,y, r,g,b)
    end
    progress (y/height)
  end
end

kuwahara(2)